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ABSTRACT 

Motivation: Several types of studies, including genome-wide associ- 
ation studies and RNA interference screens, strive to link genes to 
diseases. Although these approaches have had some success, gen- 
etic variants are often only present in a small subset of the population, 
and screens are noisy with low overlap between experiments in dif- 
ferent labs. Neither provides a mechanistic model explaining how 
identified genes impact the disease of interest or the dynamics of 
the pathways those genes regulate. Such mechanistic models could 
be used to accurately predict downstream effects of knocking down 
pathway members and allow comprehensive exploration of the effects 
of targeting pairs or higher-order combinations of genes. 
Results: We developed methods to model the activation of signaling 
and dynamic regulatory networks involved in disease progression. Our 
model, SDREM, integrates static and time series data to link proteins 
and the pathways they regulate in these networks. SDREM uses prior 
information about proteins' likelihood of involvement in a disease 
(e.g. from screens) to improve the quality of the predicted signaling 
pathways. We used our algorithms to study the human immune re- 
sponse to H1N1 influenza infection. The resulting networks correctly 
identified many of the known pathways and transcriptional regulators 
of this disease. Furthermore, they accurately predict RNA interference 
effects and can be used to infer genetic interactions, greatly improving 
over other methods suggested for this task. Applying our method to 
the more pathogenic H5N1 influenza allowed us to identify several 
strain-specific targets of this infection. 

Availability: SDREM is available from http://sb.cs.cmu.edu/sdrem 
Contact: zivbj@cs.cmu.edu 

Supplementary information: Supplementary data are available 
at Bioinformatics online. 

1 INTRODUCTION 

A wide variety of experimental and computational approaches 
have been used over the past few years to screen for genes that 
play important roles in human disease. These include RNA inter- 
ference (RNAi) screens (Mohr et al, 2010), which knock down 
human genes to quantify the phenotypic effects in diseases such 
as HIV (Brass et al, 2008) and influenza (Brass et al, 2009; 
Karlas et al, 2010; Konig et al, 2010) infection, and genome- 
wide association studies (GWAS), a powerful approach for un- 
covering strong connections between genetic variations such as 
single-nucleotide polymorphisms (SNPs) and disease traits 
(Altshuler et al, 2008). 

Although these and similar sequencing-based methods have 
been widely applied and were able to identify many relevant 
genes, they also suffer from important drawbacks that limit 
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their effectiveness. RNAi screen hits are often not reproducible 
among different labs and experimental settings. For example, the 
overlap among hits in HIV screens is low. Three genome-wide 
screens identified 842 genes that impact HIV replication, but 
only three genes were common to all screens (Bushman et al, 
2009). Similar low overlap is observed for other diseases as well 
(Stertz and Shaw, 201 1). Interpretability of screen hits remains a 
challenge, as the screens cannot explain why a gene impacts a 
disease. Similarly, interpretation of GWAS hits is not straight- 
forward as they oftentimes lie outside of coding regions of the 
genome (Schaub et al, 2012). Even those variants lying within a 
coding region only explain a small fraction of affected individ- 
uals for several diseases and conditions (Maher, 2008). 

These observations motivated methods that integrate GWAS 
or perturbation information with functional annotations and 
other types of genomic data (most notably gene expression and 
protein interaction data). For example, pathway-based GWAS 
(Wang et al, 2010) maps SNP P-values to proximal genes and 
links pathways to diseases by testing whether any predefined 
pathways are enriched for genes that flank significant SNPs. 
However, this and similar strategies are dependent on the quality 
of the annotated pathways preventing their use for important 
diseases that are poorly represented in public pathway databases 
(Wang et al, 2010). 

A variety of gene prioritization algorithms, which have been 
extensively reviewed (Moreau and Tranchevent, 2012; Piro and 
Di Cunto, 2012) and benchmarked (Bornigen et al, 2012), use 
known disease genes from literature curation, GWAS or other 
sources as seeds to search for similar genes that are likely to be 
related to the same disease. The evidence used to define gene 
similarity includes text mining, pathway membership, functional 
annotations, protein properties, sequence, co-expression and 
proximity in protein-protein interaction (PPI) networks. 
However, these methods do not attempt to mechanistically 
model the pathways that are activated during disease progression 
and response and thus do not provide a model to explain the 
various observations. Instead, they use a 'black box' approach; 
they take seed genes as input and produce a ranked list of can- 
didate genes without explaining the predicted relationship to the 
disease, which may impair further analysis (Moreau and 
Tranchevent, 2012). 

Other algorithms use perturbation and genetic data to link 
proteins to diseases, while at the same time using other genomic 
datasets and interaction networks to suggest how these proteins 
may be involved in disease response (Huang and Fraenkel, 2009; 
Kim et al, 201 1; Ourfali et al, 2007; Tuncbag et al, 2013; Yeang 
et al, 2004; Yeger-Lotem et al, 2009) However, these methods 
have only used static models so far. In such models, the pheno- 
typic outcomes of the network perturbations (or genetic variants) 
are typically differentially expressed genes. The algorithms then 
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connect possible disease sources (for example, genetic variants, 
perturbed genes, receptors or proteins that directly interact with 
a pathogen) with the differentially expressed genes and use the 
resulting source-target pathways to explain the role various pro- 
teins play in the disease response. 

Static models ignore the dynamic nature of responses and the 
temporal changes in active genes and pathways. Indeed, many 
transcriptional studies of disease highlight the importance of the 
temporal aspects, for example, the rapid immune response follow- 
ing infection (Li et al., 201 1; Shapira et al., 2009) and longitudinal 
clinical studies (Desai et al., 2011). Static modeling of dynamic 
interactions aggregates several different pathways at once, which 
leads to missing key factors that are only active at specific times 
during response (Ernst et al., 2007). In addition, many static al- 
gorithms do not attempt to infer the directionality of the edges on 
the pathways despite the importance of edge orientation in cor- 
rectly representing signaling paths (Gitter et al., 2011), 

To address these issues, we developed the Signaling and 
Dynamic Regulatory Events Miner (SDREM) (Gitter et al., 
2013). SDREM combines an input-output hidden Markov 
model (IOHMM) (Bengio and Frasconi, 1995) with a combina- 
torial optimization algorithm to reconstruct dynamic models of 
the signaling and regulatory networks using high-throughput 
data including PPI, transcription factor (TF) binding and time 
series gene expression data. The signaling pathways inferred by 
SDREM are directed cascades that connect sources and target 
TFs. These TFs are in turn used as part of the reconstructed 
dynamic regulatory network to explain the temporal gene expres- 
sion. SDREM was developed and tested on yeast and was shown 
to accurately reconstruct yeast stress response networks leading 
to new insights into the specific pathways controlling such re- 
sponses (Gitter et al., 2013). 

Although SDREM was successfully applied to model organ- 
isms, the application to human disease remained a challenge. The 
larger scale of the human interaction network makes search and 
inference much more complicated. In addition, several data 
sources available for studying diseases (including RNAi screens 
and GWAS data) were not used by the original SDREM algo- 
rithm. Finally, human networks are more complex and involve 
extensive redundant parallel pathways (Logue and Morrison, 
2012), which means that new methods are required to assess 
the ability of proteins and combinations of proteins to affect 
disease progression. 

Here we present new computational methods that extend 
SDREM by designing a new target function that uses prior 
node (protein) information, developing new optimization meth- 
ods for searching large networks and creating new gene ranking 
metrics based on target connectivity as a post-processing step. 
Combined, these extensions allowed us to apply SDREM to 
study human response to influenza A viral infection. The result- 
ing networks identified several key proteins and pathways known 
to be involved in H1N1 influenza response and predicted novel 
targets as well. To comprehensively test our method's ability to 
identify potential drug targets, we used the reconstructed net- 
works to predict the results of RNAi screens demonstrating 
that SDREM can make accurate predictions regarding such po- 
tential targets. SDREM correctly predicts more RNAi screen 
hits than top genome-wide gene prioritization algorithms, and 
its predictions are better suited for experimental validation 



because they are placed in the context of influenza-specific path- 
ways. We also used SDREM to predict genetic interactions that 
can lead to reduced viral load. Finally, we used SDREM to study 
the more lethal H5N1 influenza. Comparing the networks con- 
structed for H1N1 and H5N1 allowed us to predict genes that 
are uniquely involved in H5N1 response, shedding light on the 
specific characteristics of this infection. 

2 METHODS 

2.1 SDREM: Signaling and Dynamic Regulatory Events 
Miner 

SDREM (Gitter et at., 2013) links the two core components of the cellular 
response to an external infection or treatment: signaling cascades and 
transcriptional regulation. To reconstruct such models, SDREM starts 
with the upstream proteins that detect the invading pathogen and deter- 
mines signaling pathways that allow such proteins to communicate with 
downstream TFs. The activated TFs coordinate changes in the expression 
of their bound genes, which allows the cell to adapt to the new condition. 
SDREM integrates general TF binding data and PPI network informa- 
tion with condition-specific time series gene expression data to build such 
models. 

To find TFs that are end points of the signaling pathways initiated by 
the upstream proteins and responsible for the observed temporal tran- 
scriptional changes SDREM iteratively combines two computational 
modules. For the regulatory part of the model, SDREM uses an 
IOHMM to analyze time series expression profiles and identify split 
events — places in the time series where a group of co-expressed genes 
starts to diverge (Ernst et al., 2007). These split events are annotated 
with TFs that are predicted to control them allowing the method to 
assign temporal information to the (static) TF-gene interaction data. 

Using the putative active TFs predicted for various splits, SDREM 
assesses the likelihood that each of these TFs is indeed responding to the 
infection or treatment. A TF that is not well-connected to the upstream 
source proteins in the signaling network is unlikely to be a major driver of 
the transcriptional response. To search for such connections, SDREM 
enumerates all possible depth-bounded paths from the sources to the TFs 
in the undirected PPI network. Next, SDREM orients the PPI network 
forcing information to flow along each edge in a single direction (Gitter 
et at., 2011). To find the optimal orientation SDREM maximizes the 
following combinatorial function: 

£>0>)w(p) (1) 

where P is the set of all unique depth-bounded paths between sources and 
TFs, I s (p) is an indicator function that has the value 1 if path p is satisfied 
and w(p) is the path weight. A satisfied path is a path in which all of the 
edges have been oriented such that the TF is reachable from the source. 
Our initial edge orientation analysis (Gitter et at., 2011) used only edge 
confidence to weight paths, which in practice reduced w(p) to 

w(p) = Yi <e) (2) 

where p is a source-target path and e is an edge on the path. Edge weights 
w(e) depend on our confidence that a specific interaction exists 
(Supplementary Table SI) so that more confident pathways have 
higher weights (Gitter et al., 2011). 

SDREM uses the source-TF pathways in the oriented network to 
quantify how well the putative active TFs are connected to the sources. 
This information is fed back into the temporal gene expression analysis 
and SDREM's two stages (IOHMM and network orientation) iterate for 
a fixed number of rounds. In all rounds (except for the first when path- 
ways have not yet been learned), the IOHMM incorporates a prior that 
represents the likelihood each TF is activated by the signaling pathways 
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predicted in the previous round. In addition, the pathways can be used to 
identify other proteins (not TFs) that are involved in many high-weight 
satisfied paths and therefore important to the response. 

As mentioned above, we previously used SDREM to study yeast stress 
response. Conceptually, SDREM is a general algorithm and can be read- 
ily applied to other organisms as well. However, when attempting to use it 
to model human disease response we faced several new computational 
challenges. First, the complexity of the human interaction networks 
means that SDREM's search algorithm would be too time-consuming 
and thus required specialized approaches to improve the runtime and 
accuracy. Second, while in yeast we relied only on edge and target con- 
fidence to identify high-scoring pathways, human data provides informa- 
tion about the nodes (proteins) as well leading to a new objective function 
for finding well-connected TFs. Finally, identifying key signaling proteins 
by determining their involvement in high-scoring pathways, as was done 
for yeast, may be less important for the human network. Instead, to 
identify potential drug targets, we would like to develop methods that 
can automatically determine the effects of single or double knockdowns 
on the ability of human cells to respond to infection or treatment. Below 
we discuss the new computational methods we developed to address these 
issues. 

2.2 Algorithm parallelization and source-target pathway 
approximations 

We explored several avenues for speeding up the computation required to 
reconstruct the human signaling networks. Although it is challenging to 
parallelize the orientation step (because local changes in edge direction 
may affect other edges in different parts of the network) other aspects of 
the algorithm could be parallelized including the randomization tests to 
assess TF connectivity scores and the initial depth-first search. 

However, these and other precomputing steps (Supplementary 
Information) do not reduce the time it takes to orient the network. The 
Maximum Edge Orientation problem is NP-hard (Gitter et al., 201 1) and 
we already solve it using a heuristic approach, suggesting that a further 
approximation may not negatively impact our results. We thus modified 
the parallel path enumeration algorithm to only store the top m paths 
from any source to any TF in our dataset, ranking the paths by path 
weight. Considering only the top paths also enables us to include early 
termination in the depth-first search's branch traversal, further reducing 
runtime. Evaluating the objective function requires summing the weights 
of all satisfied paths, and for every potential edge flip that is considered 
during greedy local search we must determine which paths are still satis- 
fied. We now approximate the calculation of the orientation objective 
function by only examining these top m undirected paths. In test 
cases with millions of potential paths, the correlation between the node 
scores obtained using the exact objective and those obtained with the 
approximated objective when only using m = 100 000 was >0.999 
(Supplementary Figs SI and S2). Therefore, we set m = 100 000 for 
our H1N1 and H5N1 analysis. 

2.3 Incorporating RNAi screens 

When modeling human response, we can sometimes use additional 
sources of information regarding the involvement of a specific protein. 
Although in the original SDREM formulation (Gitter et al, 2013) all 
proteins were assumed to have the same likelihood of participating in the 
response pathways, information such as knockdown phenotypic effects or 
GWAS data can increase our prior belief that a protein participates in 
one of the response pathways. To capture such information, we can 
modify the target function discussed above (Equation 2) to incorporate 
node priors as follows: 

w(p) = w(t)Y\w(v)Y\w(e) (3) 

yep i-sp 



where p is a source-target path, t is the target on that path, v is a vertex 
on the path, e is an edge on the path and the function h>(*) is the edge 
confidence or node prior. Equation 3 attempts to find paths that contain 
proteins that are likely involved in the response based on the screen data 
as well as highly reliable protein interactions. Because the optimization 
function in Equation 2 is NP-hard (Gitter et al., 2011), the target 
function in Equation 3 is also NP-hard and the heuristic we used to 
solve the original (edge only) function can be applied to this formulation 
as well. 

Using prior information regarding protein involvement also helps us 
better resolve edge orientation in the signaling networks. Owing to the 
larger PPI and regulatory networks for human data, there are many more 
possible connections from sources to targets and disagreements about the 
orientation of individual PPI in human models. A gene that is associated 
with a phenotypic change in an RNAi screen leads to higher weights for 
the (directed) pathways in which it is contained; thus, these are preferred 
during the network orientation. Due to pathway redundancy, the con- 
verse is not necessarily true. Genes that are negative screen hits may still 
be highly relevant to signaling pathways. 

Although related approaches use the screen hits directly as sources in 
the network (Yeger-Lotem et al., 2009), we place less trust in the RNAi 
data. Independent RNAi screens can exhibit low overlap (Bushman et al, 
2009; Stertz and Shaw, 2011) in part owing to the impact of differences in 
methodology or cell population context (Snijder et al., 2012). Supplemen- 
tary Table S2 demonstrates this disagreement for RNAi screens for 
H1N1 influenza infection. No genes are hits in all five screens, and 
only a single gene is detected in four of the five screens [note that two 
of the screens (Bortz et al., 2011; Shapira et al., 2009) are targeted, not 
genome-wide]. 

Thus, to derive a prior based on screen results, we convert the RNAi 
data into vertex weights as follows: 

[ 1 - (1 -c)", if«>0 
w(y) = (4) 
I 0.5, otherwise 

where w(v) is the weight assigned to a vertex (gene), c is the confidence in 
the screen in the range [0, 1] and n is the number of screens that report v 
as a hit. We set c = 0.75 in all analyses here but could incorporate bio- 
logical knowledge to set different confidence levels for different screens. 
These node priors are used directly in the formula for path weights [u(v) 
in Equation 3] so that paths containing many screen hits have higher 
weights. 



2.4 Predicting RNAi effects 

To predict RNAi screen hits for new conditions, we developed methods 
to estimate the in silico effects of removing a protein from the signaling 
network component of an SDREM model. Instead of directly linking the 
sources and differentially expressed genes, we compute how the connect- 
ivity to the TFs is affected when a node is removed. This allows us to 
leverage the fact that each key TF affects many genes (often hundreds) so 
blocking access to such TFs significantly impacts the cell's ability to 
mount an effective response. We devised several scoring metrics to 
quantify the effect of node deletion on the targets. These metrics vary 
along three lines: (i) All versus Top denotes whether all satisfied paths or 
only the top-ranked paths are used to calculate change in connectivity, 
(ii) Source-target Pairs versus Targets determines whether a target's 
connectivity is evaluated separately for every source (i.e. each source ac- 
tivates a target differently) or if a target is considered to be disconnected 
only when it is no longer reachable from any source (i.e. any source 
can activate the target), (iii) Weighted versus Unweighted specifies if 
connectivity is treated in a binary (connected or disconnected) or con- 
tinuous (fraction of connectivity remaining) fashion. The score for the 
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weighted variant is 



E 



score w (A) -■ 



e m 
m 



(5) 



where A is the deleted node, T is the set of all targets, P(t) is the set of 
paths to the target / to be considered (depending on the choice of All 
versus Top and Pairs versus Targets), w(p) is the path weight, /(*) is an 
indicator function and N(j>) is the set of nodes on the path. The Pairs 
option replaces the summation over targets by a double summation over 
sources and targets and updates the denominator accordingly. Intuitively, 
this score is the fraction of path weight that exists along paths that can 
still activate t after A is deleted averaged across all targets. The un- 
weighted score, which reflects the fraction of targets that are still reach- 
able after removing node A, is 



score u (A) ■■ 
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(6) 



2.5 Predicting genetic interactions 

Although screen hits for individual genes are often noisy, higher-order 
knockdowns (of two or more genes) may prove to be more robust 
because they can target several pathways simultaneously. However, ex- 
perimentally testing all possible combinations, even for pairwise inter- 
actions, is not feasible given the large number of human genes. Our 
method provides a viable alternative: rather than performing an experi- 
mental screen, perform in silico knockdowns of all pairwise combinations, 
score each pair and then only experimentally test the top-ranking pairs. 

To derive a score for pairwise knockdowns, we mimicked experimental 
studies of genetic interactions in yeast (Collins et al., 2010; Jonikas et al., 
2009; Tong et al., 2004). Initial work classified gene pairs as either inter- 
acting or not with less emphasis on the strength of the interaction (Tong 
et al., 2004). More recent work has focused on quantifying genetic inter- 
actions on a continuous scale. 
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where e^s is the interaction between genes A and B and P AB is the 
phenotype when both A and B are deleted (Collins et at, 2010). 
Typically the expected phenotype is defined as the product of the pheno- 
types observed in the individual single deletions 



pobscrvcd pobserved 
r A B 



(8) 



Using this definition, negative interactions, the type we are primarily 
interested in because they indicate the pair blocks disease-related pheno- 
types, occur when the double knockdown has a stronger effect than ex- 
pected because stronger effects correspond to lower values of p^ 1 ™ 1 . 

In yeast experimental work, colony size is a typical phenotype (Collins 
et at, 2010; Tong et al, 2004) because it approximates growth rate, but 
other possibilities exist (Jonikas et a!., 2009). Single human knockdown 
screens often use viral load as the phenotype (Konig et al., 2010). In our 
simulations, the score defined in Equation 5 is the in silico phenotype 
^observed similar to colony size and viral load, in our scoring function, 
more significant deletions result in lower values, and it is meaningful to 
take the product of the scores from two individual deletions. To calculate 
poteerved^ we g enera ij zec [ Equation 5 to consider the simultaneous removal 
of two nodes from the signaling network. 
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Equation 9 represents the average fraction of path weight that remains 
after removing paths that contain node A or node B. 

2.6 Data 

The interaction network contained PPI from BioGRID (Stark et al., 
2006), post-translational modifications and PPI from the Human 
Protein Reference Database (Mishra et al., 2006) and TF-gene binding 
predictions (Ernst et al., 2010) that were processed as described in (Schulz 
et al., 2012) (Supplementary Information). The H1N1 data consisted of 
temporal gene expression (10 time points) (Shapira et al., 2009), five 
RNAi screens (Bortz et al., 2011; Brass et al., 2009; Karlas et al., 2010; 
Konig et al., 2010; Shapira et al., 2009) and 204 source proteins that 
detect influenza infection or directly interact with influenza proteins 
(Supplementary Table S3). Similarly, for H5N1, we collected temporal 
expression data (six time points) (Li et al., 2011), one small-scale RNAi 
screen (Bortz el al., 2011) and 41 sources (Supplementary Table S4). 
Source proteins were compiled from the VirHostNet database (Navratil 
et al., 2009), large-scale host-pathogen PPI studies (Shapira et al., 2009; 
Tafforeau et al., 2011) and the literature (Supplementary Information). 

We set all of SDREM's parameters to their default values (Gitter et al., 
2013) except the number of top-ranked paths used to score nodes and 
targets, which we set to 1000. 

2.7 Gene prioritization algorithms 

To benchmark SDREM's ability to predict H1N1 RNAi screen effects, 
we tested two genome-wide gene prioritization algorithms. Endeavour 
(Aerts et al., 2006; Tranchevent et al., 2008) and Pinta (Nitsch et al., 
2009, 2011). Endeavour ranks genes using many lines of evidence such 
as functional annotations, gene expression, interaction networks, text 
mining and sequence similarity and combines the individual rankings to 
create a global prioritization. Pinta is designed for diseases where a set of 
seed genes is not available but there are differential gene expression data 
between healthy and diseased individuals. It ranks genes based on the 
expression levels of their neighbors in an interaction network. The human 
proteins that interact with H1N1 proteins and differentially expressed 
genes were provided as input to both Endeavour and Pinta. We used 
EDGE (Leek et al., 2006) to identify significantly differentially expressed 
genes by comparing the mock treatment and viral treatment time courses 
and using 500 null iterations. Endeavour does not use weights on the 
input genes so we provided the sources and all genes that had the most 
significant q-value from EDGE (1.64E-6). We performed a genome-wide 
ranking using all lines of evidence. For Pinta the weight of a gene was 
— log 10 (<7— value), and all genes were used as input. Sources were given the 
same weight as the most significantly differentially expressed genes. All 
default settings were used. 

3 RESULTS 

Pathogens infecting cells provide a clear set of sources that ini- 
tiate the subsequent signaling and transcriptional response. In 
particular, many viruses encode only a small number of proteins 
allowing us to generate specific models for the host response that 
is triggered by host proteins that detect viral RNA or interact 
with viral proteins. 

To test SDREM's application to such responses, we focus on 
influenza A viruses because of the rich datasets available and 
their importance to global health. The 2009 swine-origin H1N1 
virus outbreak received great public attention and was declared a 
global pandemic in June of that year (Zhang el al., 2010). More 
recently, research concerning mutations of avian H5N1 influenza 
that could allow it to be transmitted among mammals via aero- 
sols has sparked immense controversy (Berns et al., 2012), 
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highlighting the threat influenza A viruses pose and the need to 
better understand their interaction with the human immune 
system. 

3.1 H1N1 influenza model 

The host response to H1N1 influenza is the best profiled in terms 
of protein interactions, functional screens and transcriptional 
effects. We leveraged these data, in particular the temporal 
gene expression, to construct a dynamic model of the human 
immune response to H1N1 infection. SDREM identified the 
TFs that control this immune response and the signaling path- 
ways that activate them — 36 internal proteins that connect 33 
target TFs to upstream nodes in the signaling network (Fig 1 
and Supplementary Table S3). These include many proteins 
known to be involved in immune response such as STAT1 
(Shuai and Liu, 2003), seven IRF family members (Honda and 
Taniguchi, 2006), three NFKB variants (Baeuerle and Henkel, 
1994), RELA (Ouaaz et al., 1999) and XBP1 (Martinon et al., 
2010). Interestingly several cancer- re la ted proteins such as AR, 
BRCA1, MYC, SMAD3, SMAD7 and TP53 appear as well. 

We used DAVID (Huang et al, 2009) to compute the Gene 
Ontology (GO) (Ashburner et al, 2000) biological process en- 
richment of the proteins predicted by SDREM, leaving out the 
sources because they are already known to be relevant to H1N1 
and would bias the results. The most significant GO terms are 
dominated by processes related to transcriptional regulation due 
to the prevalence of target TFs in our predictions. However, 
beyond these are many highly relevant enriched terms including 
'immune system development' [Benjamini-Hochberg corrected 
P-value 2.18 E-5 (Benjamini and Hochberg, 1995)], 'response to 
virus' (P= 0.0169), 'virus-host interaction' (P = 0.0289) and 
'immune response' (P= 0.0431). 

3.2 Predicting RNAi screen hits 

SDREM can rank human genes for their involvement in mediat- 
ing host response to viral infection. Even for viruses for which 
genome-wide screens are available SDREM's models are useful 
since they provide mechanistic explanations for genes' involve- 
ment in the response. In addition, the large disagreement among 
genome-wide screens (Supplementary Table S2) makes such 
models important in order to distinguish real hits from noise. 
Finally, several pathogens, including H5N1, are challenging to 
work with because they require a biosafety level 3 lab, making it 
difficult to generate genome-wide RNAi screens. Using SDREM 
to predict RNAi screen hits allows researchers to prioritize 
candidate H5N1 targets, which could later be validated 
experimentally. 

To determine which of the scoring metrics described in Section 
2.4 is most predictive of RNAi screen hits, we ran SDREM on 
the H1N1 data holding out the RNAi screen data from 
SDREM's input. We used the metrics to rank all 252 non- 
target proteins in the model using the number of high-confidence 
paths that use the node and the network degree to break ties in 
the ranking. Given the rankings for each metric, we computed 
the percentage of correct hits within the top X predictions (where 
Granges from 10 to 100). 

As seen in Table 1, SDREM performed exceptionally well on 
this task. To illustrate this, consider the number of correct screen 
hits in the top 50 ranked genes. Roughly 20 of these (depending 



on the actual metric used), or 40%, are known hits. Because all 
the 14 334 proteins in the interaction network are included in 
SDREM's search, the expected number of known hits in a ran- 
domly selected set of 50 genes is 3.1. Thus, using SDREM we 
obtain a 6.5-fold enrichment in the number of correct hits, indi- 
cating that such a method can be effectively used to design spe- 
cific experiments for other viruses as well. Similar enrichments 
are seen for other values of X. 

The best-performing metric uses only the top paths, allows 
targets to be activated by any source and uses the weighted 
score. Intuitively, including the lower confidence paths can 
hurt predictive performance because these source-target connec- 
tions may contain false-positive PPI and not actually affect con- 
nectivity when they are broken. Our results also suggest that a 
target can function if it remains connected to any source. Due to 
the many parallel paths in the signaling network, it is uncommon 
for a single node removal to completely disconnect a target. The 
weighted variant distinguishes between deletions that do not 
impact a target at all and those that remove many paths to the 
target but do not fully disconnect it. This metric's predictions 
significantly overlap with the known screen hits at all thresholds, 
which is also the case for most of the other metrics. Because the 
overlaps are significant even for the worst-performing metric 
when 20 or more genes are predicted, we conclude that the 
SDREM model itself is a powerful filter for predicting screen hits. 

3.3 Comparison with gene prioritization algorithms 

Having demonstrated that SDREM produces highly accurate 
rankings of putative RNAi screen hits, we examined whether 
existing tools could achieve similar performance. Although 
other algorithms for connecting sources and transcriptional ef- 
fects provide candidate pathways (Huang and Fraenkel, 2009; 
Kim et al., 2011; Ourfali et al., 2007; Schaefer et al, 2013; 
Tuncbag et al., 2013; Yeang et al., 2004; Yeger-Lotem et al., 
2009), they do not rank the pathway members or quantify the 
effect of knocking them down so it is impossible to compare our 
knockdown prediction results with these methods. Instead, we 
assess Endeavour (Aerts et al., 2006; Tranchevent et al., 2008) 
and Pinta (Nitsch et al, 2009, 2011), the only two gene priori- 
tization web servers in a recent benchmark (Bornigen et al., 2012) 
that can prioritize the entire genome and take seed genes instead 
of disease keywords as input. 

We ran Endeavour and Pinta using the differentially expressed 
genes identified by EDGE (Leek et al, 2006) and the H1N1 
sources as input and evaluated the number of correct screen hit 
predictions at the same thresholds used previously (Table 2). 
Because the web servers do not provide their full ranked list of 
genes, we cannot calculate AUC. SDREM outperforms the gene 
prioritization tools at all thresholds except for the top 20 
threshold where Pinta makes the same number of correct predic- 
tions and Endeavour makes two more. SDREM's advantage is 
greatest when assessing the top 100 predictions, at which point 
SDREM correctly predicts nearly twice as many RNAi screen 
hits as Pinta and 24% more than Endeavour. When using alter- 
native methods to select the input genes for Endeavour and 
Pinta, SDREM's strengths are even more pronounced 
(Supplementary Information, Supplementary Table S5). These 
results are especially impressive considering Endeavour uses 
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Fig. 1. The SDREM H1N1 response model, (a) The regulatory paths summarize the temporal patterns of the differentially expressed genes. The x-axis is 
time and the y-axis is log 2 fold change. Split events, green nodes where a regulatory path branches, are annotated with the TFs that are predicted to 
activate or repress the genes at that time point. These annotations are placed on the path immediately after the split to indicate whether the TF controls 
the upper or lower path out of the split, (b) The signaling paths from sources (red) through internal nodes (blue) to the active TFs (green). Sources 
directly interact with viral proteins or detect viral presence. The active TFs are the same TFs shown on the regulatory paths. Diamonds are RNAi screen 
hits. Solid edges are PPI whose orientation has been inferred by SDREM. Dashed edges are post-translational modifications and TF-gene binding 
interactions, which already have a known orientation 



additional input that SDREM does not including prior literature 
via text mining. These results stress the importance of building 
mechanistic models that explain why differentially expressed 
genes are affected by infection instead of directly finding genes 
that are similar to the differentially expressed genes. 



3.4 Predicting host proteins affecting H5N1 viral load 

Having successfully used SDREM to predict H1N1 screen hits 
we turned to the independent H5N1 data. Using the scores from 
the best metric, we ranked all proteins in the H5N1 SDREM 
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Table 1. The scoring metrics that were used to rank H1N1 screen hits 



Paths used 


Connectivity 


Score 


AUC 


Hits in top 10 


Hits in top 20 


Hits in top 50 


Hits in top 100 


Top 


Targets 


Weighted 


0.722 


6 (1.97E-5) 


8 (3.44 E-5) 


18 (3.24 E-9) 


42 (9.42 E-23) 


Top 


Pairs 


Weighted 


0.717 


3 (2.87 E-2) 


7 (2.88 E-4) 


20 (4.59 E-ll) 


40 (8.68 E-21) 


Top 


Pairs 


Unweighted 


0.716 


3 (2.87 E-2) 


8 (3.44 E-5) 


20 (4.59 E-ll) 


37 (5.59 E-18) 


All 


Targets 


Unweighted 


0.711 


2 (0.153) 


5 (1.09 E-2) 


20 (4.59 E-ll) 


39 (7.82 E-20) 


All 


Targets 


Weighted 


0.706 


3 (2.87 E-2) 


6 (1.97E-3) 


18 (3.24 E-9) 


39 (7.82 E-20) 


Top 


Targets 


Unweighted 


0.704 


2 (0.153) 


6 (1.97E-3) 


19 (4.02 E-10) 


39 (7.82 E-20) 


All 


Pairs 


Weighted 


0.702 


3 (2.87 E-2) 


6 (1.97E-3) 


18 (3.24 E-9) 


36 (4.43 E-17) 


All 


Pairs 


Unweighted 


0.676 


2 (0.153) 


6 (1.97E-3) 


18 (3.24 E-9) 


36 (4.43 E-17) 



Note: The metrics are sorted by area under the curve (AUC). The number of known screen hits recovered at various thresholds is shown with the significance (in parentheses) 
calculated using Fisher's exact test. 



Table 2. Comparison of SDREM, Endeavour and Pinta gene rankings 



Algorithm 


Settings 


Hits in 
top 10 


Hits in 
top 20 


Hits in 
top 50 


Hits in 
top 100 


SDREM 


Top, targets, weighted 


6 


8 


18 


42 


Endeavour 


All evidence 


5 


10 


17 


34 


Pinta 


Default 


5 


8 


12 


22 



model (Supplementary Table S4) and compared the ranks with 
those we obtained using the same scoring metric on the H1N1 
data. We also examined the degree of the top-ranked predictions 
because we expected that high-scoring nodes would have high 
degree because such nodes are likely to affect a large number of 
targets when deleted. Table 3 lists the top 25 H5N1 predictions. 
Recall from Equation 5 that the scores denote the fraction of 
target connectivity remaining after the in silico deletion so lower 
scores translate to a stronger effect. Six of the H5N1 predic- 
tions— STAT3, CASP8, HSF1, ERBB3, NRIP1 and 
PRMT1 — are particularly interesting because they are neither 
known H1N1 RNAi screen hits nor in the top 100 H1N1 pre- 
dicted hits. Two of these, CASP8 and ERBB3, have been re- 
ported to directly interact with H5N1 viral proteins but not 
H1N1 proteins even though the HINl-human PPI have much 
greater coverage, suggesting that they may indeed play distinct 
roles. Several top-ranked H5N1 proteins are sources (directly 
interacting with the viral proteins) or high-degree nodes. In con- 
trast, NRIP1 and PRMT1 are neither sources nor of high degree. 
Their inclusion in the top predictions is noteworthy because the 
paths through these nodes affect targets to a greater extent than 
expected. NRIP1, also known as RIP140, is involved in inflam- 
matory response as a coactivator for NFKB (Zschiedrich et ah, 
2008). PRMT1 methylation of STAT1 is one of the ways in 
which STAT1 is regulated in the immune system (Shuai and 
Liu, 2003). 

3.5 Screen hits improve the SDREM H1N1 model 

SDREM can predict screen hits, but when RNAi data are al- 
ready available they do indeed lead to more accurate SDREM 
models. Here we refer to the model reconstructed using all data 
as the 'original' model and the model obtained after omitting the 



RNAi data as the 'no- RNAi' model. Of the 69 signaling proteins 
and TFs predicted in the original SDREM H1N1 model, 38 
(55%) are screen hits versus only 20 of 61 predictions (33%) in 
the no-RNAi model. The original and no-RNAi models overlap 
substantially. However, some important immune responders, 
including IRF2 and NFKB1, were only detected in the original 
SDREM model. 

The original SDREM model better corresponds to related 
annotated pathways as well. We used DAVID to examine the 
enrichment of Biocarta pathways (Nishimura, 2001) for SDREM 
model members. The immune-related pathways 'Human 
Cytomegalovirus and Map Kinase Pathways', 'The informa- 
tion-processing pathway at the Interferon-beta enhancer', 'T 
Cell Receptor Signaling Pathway' and 'Toll-Like Receptor 
Pathway' are enriched in the original model (all with corrected 
P<0.05) but not the no-RNAi SDREM model. 

3.6 Predicting genetic interactions 

Genetic interactions are functional interactions between pairs of 
genes where simultaneous double deletion has a smaller or 
greater than expected effect. Humans have tens of thousands 
of genes, making it impossible to comprehensively screen for 
all possible genetic interactions in a condition of interest as is 
done to identify the phenotypic effects of single gene loss. In 
contrast, the in silico analysis performed by SDREM can be ex- 
tended to pairwise or higher-order analysis. Given its perform- 
ance on the single knockdown prediction task, the ranked list of 
pairs identified by SDREM can serve as a starting point for 
follow-up validation experiments. Considering the rather disap- 
pointing agreement between screen hits performed by different 
labs, such higher-order analysis may be required to accurately 
and robustly identify targets that can impact disease progression 
and viral load. 

Based on our results for predicting H1N1 screen hits (Table 1), 
we again consider only the top-ranked paths and allow targets to 
be activated by any single source. We used our method to predict 
genetic interactions that affect H1N1 (Table 4) and H5N1 
(Supplementary Table S6) infection. Condition-specific experi- 
mental validation is necessary to confirm that these pairs form 
genetic interactions that impact viral infection, but annotated 
pathways indicate that some of these pairs do act in parallel in 
other settings. For example, RBI and TP53 are on parallel paths 
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Table 3. The top-ranked H5N1 RNAi screen hit predictions alongside H1N1 RNAi rankings and the number of screens reporting known hits 



Gene 


H1N1 source 


H 5N1 source 


Degree 


H1N1 RNAi 


H5N1 RNAi 


H5N1 score 


H1N1 rank 
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HSF1 


N 


N 


217 


0 


0 


0.950 


N/A 


16 


EP300 


N 


N 


377 


1 


0 


0.951 


3 


17 


BRCA1 


N 


N 


301 


0 


0 


0.954 


49 


18 


NUP98 


N 


Y 


36 


2 


0 


0.955 


N/A 


19 


ERBB3 


N 


Y 


37 


0 


0 


0.963 


N/A 


20 


NRIP1 


N 


N 


48 


0 


0 


0.964 


N/A 


21 


STAT1 


N 


N 


642 


0 


0 


0.964 


22 


22 


PRMT1 


N 


N 


70 


0 


0 


0.964 


147 


23 


KPNA1 


Y 


Y 


26 


1 


1 


0.967 


216 


24 


HSP90AA1 


Y 


N 


144 


2 


1 


0.968 


9 


25 



Note: N/A indicates that the gene was not included in the SDREM H1N1 model. 



Table 4. The top 10 predicted H1N1 genetic interactions 



Gene A 


Gene B 


(AB 


pob 
r AB 


pex 
r AB 


pob 


pob 
' B 


EP300 


TP53 


-0.0077 


0.8152 


0.8229 


0.9158 


0.8986 


TRAF2 


UBE2I 


-0.0070 


0.8275 


0.8345 


0.9348 


0.8927 


UBC 


UBE2I 


-0.0070 


0.8256 


0.8326 


0.9327 


0.8927 


RBI 


TP53 


-0.0068 


0.8316 


0.8384 


0.9330 


0.8986 


TP53 


TRAF2 


-0.0066 


0.8333 


0.8400 


0.8986 


0.9348 


RBI 


UBE2I 


-0.0057 


0.8272 


0.8329 


0.9330 


0.8927 


EP300 


UBC 


-0.0057 


0.8485 


0.8541 


0.9158 


0.9327 


EP300 


TRAF2 


-0.0055 


0.8506 


0.8561 


0.9158 


0.9348 


EIF2AK2 


UBE2I 


-0.0053 


0.8432 


0.8485 


0.9505 


0.8927 


NPM1 


UBE2I 


-0.0052 


0.8442 


0.8494 


0.9515 


0.8927 



in the Biocarta pathway 'Tumor Suppressor Arf Inhibits 
Ribosomal Biogenesis' as are TP53 and TRAF2 in KEGG's 
(Kanehisa and Goto, 2000) 'Pathways in cancer'. Other predicted 
pairs such as ILF3 and PA2G4 (fifth on the H5N1 list) are es- 
pecially interesting because the two proteins are not high-degree 
nodes in the PPI network. 

4 DISCUSSION 

SDREM is unique in that it combines time series and static data 
to model both signaling and dynamic regulatory networks. This 



allows it to infer, more confidently, the TFs that are at the end 
points of signaling cascades and control human disease response. 
To manage the complexity of human interaction networks, we 
extended SDREM to leverage gene priors from RNAi screens in 
its objective function and incorporated several algorithmic im- 
provements. As we have shown for influenza infection, by recon- 
structing disease response networks we can accurately identify 
key signaling pathways and nodes (proteins). SDREM can be 
applied successfully in many different settings even when the PPI 
data, node priors or gene expression time points are sparse 
(Supplementary Information, Supplementary Tables S7 and S8). 
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Given the predicted directed signaling pathways, we can esti- 
mate the phenotypic effects of knocking down genes by assessing 
how strongly the downstream TFs are affected. These techniques 
can be used to infer putative drug targets for hard-to-study con- 
ditions (such as H5N1 infection) and for combination of targets, 
which can lead to more robust treatments. Successfully predict- 
ing the effects of pairwise and higher-order gene knockdowns 
can guide targeted experimental validation, a major contribution 
because exhaustive pairwise RNAi screening is currently infeas- 
ible and signaling pathway redundancy can limit the effectiveness 
of drugs that target individual genes (Logue and Morrison, 
2012). Owing to this pathway redundancy and false negatives 
in the existing RNAi screens, the precision we calculate for 
SDREM's RNAi effect predictions (Table 1) is conservative. 
Many proteins predicted by SDREM (for example STAT1, 
RELA, NFKB2 and several IRF TFs) are not screen hits but 
are known to be important to the immune response. 

An important advantage of SDREM over previous RNAi 
screen and host-pathogen PPI studies is the ability to infer 
new pathways from general interaction data. Previous studies 
[for example, (Shapira et al., 2009)] and pathway-based GWAS 
(Wang et al., 2010) often rely on known curated pathways that 
are incomplete and not always relevant to the disease studied. In 
contrast, using condition-specific time series data and sources, 
SDREM can predict which candidate proteins are involved in 
the signaling pathways and which are not. 

Unlike the gene scores from Endeavour, Pinta and other gene 
prioritization algorithms that are only defined for single genes, 
SDREM predicts functional disease genes by modeling the 
impact of removing a gene from the disease-specific signaling 
pathways. This simulated phenotype can be naturally extended 
to pairs of genes, allowing SDREM to predict genetic inter- 
actions analogously to experimental approaches. Existing genetic 
interaction prediction algorithms require a partial set of known 
genetic interactions (Bandyopadhyay et al., 2012; Qi et al., 2008; 
Wong et al., 2004; Zhong and Sternberg, 2006), which prevents 
their application in human. Furthermore, SDREM's genetic 
interaction predictions are condition-specific. 

Our influenza analysis focused on data from human cell lines, 
but some of SDREM's most exciting future applications will 
involve data from individual patients. We would like to use the 
extensions presented in this article to enable GWAS data to be 
used as node priors. One possible approach we intend to explore 
is to map SNP P-values to gene scores as in pathway-based 
GWAS (Wang et al., 2010). In addition, SNPs in non-coding 
regions with potential regulatory functions could be used to sug- 
gest which TFs' binding is disrupted (Schaub et al., 2012), pro- 
viding priors for both the network orientation and the temporal 
expression analysis. 
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